/* visFlowTopology.h
*/

#ifndef VIS_FLOW_TOPO
#define VIS_FLOW_TOPO

#include "visStreamline.h"
#include "vis2by2RealMatrixEE.h"

namespace visualization
{
	// The data field processed here is two dimensional.
	static const int dim = 2;

// The struct is used to save the information about saddles.
struct saddlePoint
{
	/* The coordinates of a saddle.
	 * This is added after the function "groupSaddles" is added, because the function
	 * can output a saddle point which is not located just at the position of a vertex,
	 * so using id will be impossible to determine the location of a saddle. Therefor,
	 * the variables x, y and z are added here.
	 * What's more, the variable id will be commented and not used any more.
	 *                                  2008 - 9 - 16. by lu.
	 */
	double x, y, z;

	/* A 2 by 2 array to save eigenvectors.
	 * [v11 v12 v21 v22]
	*/
	double eigenvec[dim * dim];
};

enum Type{ATTRACTINFOCUS, REPELLINGFOCUS, ATTRACTINGNODE, REPELLINGNODE, CENTER};
struct criticalPoint
{
	Point position;
	Type type;
};

class visFlowTopology
{
public:
	visFlowTopology(Model *m) : dataField(m), no_zero_point(1)
	{
		integrationForwardCurves.setDataField(dataField);
		integrationBackwardCurves.setDataField(dataField);
	}
	visFlowTopology() : no_zero_point(1){}
	void setDataField(Model *m)
	{
		dataField = m;
		integrationForwardCurves.setDataField(dataField);
		integrationBackwardCurves.setDataField(dataField);
	}
	~visFlowTopology()
	{
	//	delete []ForwardSources;	
	//	delete []BackwardSources;	
	}

	/* True is returned if the input vector is zero and in the parameter "pos" 
	 * the position of the zero point is returned. Otherwise, flase is returned.
	 */
	bool isZeroPoint(size_t id, double *pos);

	/* 0 is returned if the input jac is a zero matrix, otherwise
     * 1 is returned.
	*/
	bool matrixIsZero(const double jac[dim * dim]);

	/* Compute the Jacobian matrix at the point located at pos. If succussfuly computed,
	 * true is returned, otherwise false is returned.
	*/
	bool computeJacobian(double *pos);

/*	// Calculate the eigenvalues and eigenvectors.
	void calEigenvalAndVec();
*/
	// Classify the zero points according to their eigenvalues.
	bool classify(const double *pos);
	
	// This is for saddle.eigenvector only.
	void normalizeVec(double vec[4]);

	/* This function will search for zero points and save them to corresponding 
	 * variables such as repellFocus, attractingFocus. 
	*/
	//1. For each point in dataField
	//2.    if(p(v)==0)
	//3.      {
	//4.		no_zero_point = 0;
	//5.		Jacobian;
	//6.		eigenvalue: point type;
	//7.		eigenvector;
	//8.       }
	void determineCriticalPoint();

	/* For each saddle, select integration starting point at a specified distance away from the 
	 * saddle along its eigenvectors' directions.
	 *
   * The integration curves extends until they reach another critical point or the boundry.
	 * Both of the two cases can be determined automatically by the visRungeKutta4 integration 
	 * method implemented before, so what's needed is just to integrate from the saddles without 
	 * extral concerns.
	*/
	void integrateCurves(intgResults *&forward, intgResults *&backward);

	// Calculate the eigenvalue and eigenvector respectively. Note that the
	// dimension of the jac should be given, and the returned results are
	// saved in the array eigval and eigvec.
//	bool calculator(double jac[], complex_ eigval[], complex_ eigvec[]);

	// Get the sources.
//	visualization::Point *getSources() {	return sources;	}

	void setRdAndSmpRate(double rds, double smp) {	radius = rds;	smpRate = smp; }
	void setZero(double z = 1e-5) { zero = z; }
	void setDelt(double d = 1e-2)	{	delt = d;	}
	void setDiavationFromSaddle(double d) {	diavation = d; }
	void setNeighbourNum(int n = 10) {	neighbourNum = n;	}
	void setDisError(double de = 1) {	disError = de; }


	/**
	 *
	 *
	 */
	void preProcessSaddles();

	vector<saddlePoint> getSaddles(size_t &num) { num = saddles.size();	return saddles; }

	void getSources(Point *&forward, Point *&backward)
	{
		forward = ForwardSources;
		backward = BackwardSources;
	}

/**
 * The variable is set as protected to facilitate it's usage in the class visTopoClustering
 */
protected:
	// The datafield.
	Model *dataField;
	vector<saddlePoint> saddles;
	vector<criticalPoint> criticalPoints;

	Point *ForwardSources;	
	Point *BackwardSources;	
	
  /** As a result of discretization of the problem, one saddle point in the data field may lead to
	 * many ones that can be detected by the program.So, it's necessary to group the saddles that actually 
	 * belong to the same one saddle point.
	 *
	 * This function clusters the saddle points saved in the variable saddles, calculates new locations and
	 * of eigenvectors of saddle points that turn out to be the same one saddle according to some metric.
	 *
	 * There is one solution to the combination:
	 * 1. find the two saddle point that are sufficient close to each other;
	 * 2. compare their eigenvectors and the ones with differences of eigenvectors under some threshold
	 *    is considered as simply the same saddle point.
	 *
	 * To group the saddles, calculate the distance between them, and two points that are enough 
 	 * close to each other and have almost the same eigenvectors are considers as one saddle point.
   * The distance can be two times radius which is used to determine zero points.
	 */
	void groupSaddles(void);

	void removePesudo(void);

	void adjustSaddlePosition(void);
	
	/**
	 * Shift the 'pos' by dist along each direction of the eigenvectors of saddle 'sdl', and 
	 * compute and return the normalized projections.
	 */
	void projection(const double pos[3], const size_t sdl, const double dist, double prj[4]);

private:
	// This is used to save the integration curves.
	visStreamline integrationForwardCurves;
	visStreamline integrationBackwardCurves;

	//The matrix is stored in a 
	//colume-wise manner.
	double jacobian[dim * dim];

	// The eigenvalue and Eigenvectors of type of complex_ used by clapack.
	// The vectors are stored in a column-wise order.
	//complex_ eigenvalue[dim];
	//complex_ eigenvector[dim * dim];
	complex<double> eigenvalue[dim];
	complex<double> eigenvector[dim * dim];

	/* This is used to indicate whether there is zero point in the 
	 * data field.
	*/
	bool no_zero_point;

	/* The paramters are used to determine whether a point is a 
	 * critical point. When doing so, a circle is made around the
	 * cadidate point, then the flow through the circle is calculated.
	 * The cadidate point may be a zero point if the resulted number
	 * is zero.
	 * @radius : the radius of the circle used.
	 * @smpRate : the number of the smpled points on the circle.
	 */
	double radius;
	double smpRate;

	void push(vector<criticalPoint> &p, const double *pos, Type t)
	{	
		Point tmp(pos[0], pos[1], pos[2]);
		criticalPoint cp = {tmp, t};
		p.push_back(cp);
	}
	void push(saddlePoint &p, const double *pos)
	{	p.x = pos[0];	p.y = pos[1];	p.z = pos[2];	}

	/* Since in the subroutine it is needed to judge whether a number is zero, so here a 
	 * uniform definition of zero is given.
	*/
	double zero;

	// When computing Jacobian matrix, a differential computation is needed. 
	double delt;// = 1e-3;

	// Parameters used for interpolation.
	int neighbourNum;
	double disError;// = 4;

	/* When computing the integration curves stemming from the saddles, it is neccessary 
	 * to begin them from somewhere apart from the saddles under consideration, since the 
	 * "velocity" at the saddles is zero. This constant is simply the distance that will 
	 * be used.
	*/ 
	double diavation;// = 0.01;

	size_t wunai;
};

}
#endif

